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Abstract 

We consider the use of a multigrid method with central differencing to solve the 
Navier-Stokes equations for hypersonic flows. The time-dependent form of the equations 
is integrated with an explicit Runge-Kutta scheme accelerated by local time stepping and 
implicit residual smoothing. Variable coefficients are developed for the implicit process 
that remove the diffusion limit on the time step, producing significant improvement in 
convergence. A numerical dissipation formulation that provides good shock-capturing 
capability for hypersonic flows is presented. This formulation is shown to be a crucial 
aspect of the multigrid method. Solutions are given for two-dimensional viscous flow 
over a NACA 0012 airfoil and three-dimensional viscous flow over a blunt biconic. 
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Introduction 

At the present time there is a strong interest in high-speed flight vehicles. Some 
examples of these vehicles are the high-speed civil transport (HSCT) and the hypersonic 
flight configurations being considered for the National Aero-Space Plane (NASP). In the 
case of the NASP, one encounters complex high Mach number phenomena and interac- 
tions, which can involve strong shock and expansion waves, in not only the flow over 
the vehicle, but also in the flow through the engines, where chemical reactions occur. 
An effective design method for such vehicles will obviously require both detailed experi- 
mental data as well as flexible, efficient, and accurate computational techniques. Robust 
prediction methods (i.e., those that can be applied on a routine basis) are not currently 
available for hypersonic flows. However, due to the significant progress during the last 
decade in the development of effective algorithms for subsonic and transonic flows, there 
are a number of opportunities for constructing improved schemes for high-speed flows. 

One powerful approach for the numerical solution of partial differential equations, 
which has been successfully applied to fluid flow problems, is multigrid. Multigrid meth- 
ods were first developed for elliptic equations. These were later extended to hyperbolic 
equations such as the time-dependent fluid dynamic equations for subsonic and tran- 
sonic flow [1, 2, 3, 4, 5]. Even for transonic cases, the steady state can retain many of 
the properties of an elliptic equation when the region of supersonic flow is limited. We 
shall show that with proper care the multigrid method still works for hypersonic flow. 
Gustafsson and Lotstedt [6] have pointed out that hyperbolic multigrid works by two 
different processes. For the long waves, the advection process is most important and 
multigrid achieves its efficiency by allowing the use of larger time steps on coarser grids. 
Hence, it is important that the smoother use large time steps. However, for the shorter 
waves, dissipation is more important and the efficiency of multigrid is based on principles 
similar to that for elliptic equations. 

Several investigators have applied multigrid to high-speed flows with varying degrees 
of success. For example, in [7] the Euler equations are solved with and without chemistry 
using the Roe scheme for spatial discretization and the ADI scheme for time marching. 
A factor of four decrease in computational time was obtained with the multigrid method 
for a simple one-dimensional nozzle flow (exit Mach number of about 3.7). A calculation 
for a Mach 22 wedge flow showed the basic scheme to be noticeably faster than the 
multigrid scheme. One notable conclusion of this work was that the performance of the 
multigrid is driven by the fluid dynamics and not the chemistry (at least for the case of 
simple reaction models). The high level of performance and widespread application of 
multigrid algorithms with central differencing and an explicit multistage time-stepping 
scheme have provided strong encouragement to make them work for hypersonic problems. 
An initial effort [8] to apply this type of algorithm resulted in numerical difficulties that 
prevented the calculation of two-dimensional flows (i.e., blunt body and wedge type) with 
a Mach number higher than about. 7. In order to compute such flows, a low Courant- 
Friedrichs-Lewy (CFL) number was required. Thus four and five stage schemes were not 
practical, since there is substantial deterioration in the high frequency damping of the 
scheme due to the large reduction in the CFL number. The CFL restriction reduced the 
potential of the scheme as a viscous flow solver. More recently an algorithm utilizing a 
semicoarsening technique, a symmetric TVD formulation, and a three stage Runge-Kutta 
scheme [9] was proposed and used to compute high Reynolds number (laminar) Mach 
10 flow over an airfoil at 10 degrees angle of attack. A good resolution of the bow shock 
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wave and a reasonable convergence rate were obtained. The method of semicoarsening 
considered required a much more complicated cycle strategy than that employed with 
standard multigrid methods. In addition, it appears to be somewhat cumbersome to 
implement in three dimensions. 

It is our contention that standard multigrid techniques can be used in conjunction 
with central differencing for hypersonic flows. To do this one must ensure that the 
basic algorithm exhibits good damping of high frequencies, on both the fine and coarse 
meshes, in the neighborhood of strong discontinuities. It becomes more difficult to 
eliminate these high frequency oscillations as the Mach number increases, since the jumps 
across shocks become larger. Thus a considerable part of the following discussion will 
concentrate on the smoothing algorithm. The fundamental features of the multigrid 
process (i.e., Full Approximation Storage scheme, grid transfer operators, fixed cycle 
strategy) are fairly standard. Other aspects, such as type of coarse grid correction 
scheme and procedure for smoothing of coarse grid corrections, found crucial in the 
present work will be emphasized. In this paper we consider a Runge-Kutta scheme [10] 
as the smoother for the multigrid method. Central differences for spatial approximations 
are augmented by an artificial viscosity based on TVD principles [ 11 ]. Several changes 
are made to the numerical algorithm so that a converged solution can be obtained for 
high-speed flows. Initially, we describe the basic Runge-Kutta method for the central- 
difference scheme with the numerical viscosity. We discuss an implicit procedure for 
accelerating the convergence of the basic scheme. This procedure removes the diffusion 
limit on the time step, which can be quite severe at hypersonic speeds. We then apply 
the multigrid method to solve both two- and three-dimensional viscous flows. 

Basic Scheme 

The basic elements of the scalar dissipation model considered in this paper were 
first introduced by Jameson, Schmidt, and Turkel [10] in conjunction with Runge-Kutta 
explicit schemes. The spatial discretization is based on central differences with an addi- 
tional artificial viscosity. This algorithm has been used by many investigators to solve the 
Euler equations numerically for a wide range of fluid dynamic applications. The same 
type of spatial discretization has been applied to alternating direction implicit (ADI) 
schemes [ 12 ] and LU factored implicit schemes [13]. In this section the basic scheme 
is briefly reviewed. Since the elements of the two-dimensional and three-dimensional 
schemes are essentially the same, we will take, for the purpose of simplicity, a two- 
dimensional point of view in describing the basic scheme. 

Consider the Euler equations in the form 

Wf + f x + g y = 0 , ( 1 ) 

where the four-component vector of conserved variables 

W = [p pu pv pE ] T , (2) 

and /, g are the corresponding flux vectors. The quantity p is the density, u and v 
are the Cartesian velocity components, and E is the specific total internal energy. The 
independent variables are time t and Cartesian coordinates (x,y). If ( 1 ) is transformed 
to arbitrary curvilinear coordinates £ = £( x,y ) and 77 = t/), then we obtain 

(J~'W) t + Ft + G n = 0 (3) 
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where J 1 is the inverse transformation Jacobian, and 


F — fVv ~ 9 X G = gx£- fy v 

In a cell-centered, finite- volume method, (1) is integrated over an elemental volume in 
the discretized computational domain, and J -1 is identified as the volume of the cell. 
Equation (3) can also be written as 

J’ l W t + AWs + BW v = 0 


where A and B are the flux Jacobian matrices defined by A -- dF/dW and B = dG/dW . 

To advance the scheme in time we use a multistage scheme. A typical step of a 
Runge-Kutta approximation to (3) is 


w (k) = IT (0) - a k At + D v G {k ~ l) - ad] , (4) 


where D ^ and D v are spatial differencing operators, and AD represents the artificial 
dissipation terms. The derivatives of the fluxes are approximated by central differences. 
Hence, we need to evaluate F at (i + j). For the physical diffusive fluxes the approxima- 
tions follow directly. Since the dependent variables are given at (z, j) we need to average 
them to calculate the convective flux at the cell face. This averaging can be done in 
several ways. In the original scheme the conservative variables were averaged to find the 
variables at the cell face and then the flux was evaluated. Moreover, in two dimensions 
five quantities were averaged, the four conservative variables and pE + p . This was done 
because enthalpy damping was used. Alternatively one can average density, momentum, 
and total enthalpy. This might be more accurate for the Euler equations since the total 
enthalpy is constant in the steady state. An alternative to this procedure is to average 
the fluxes rather than the dependent variables, where the fluxes are defined as for (3) and 
the transformation derivatives are evaluated at the cell interface. Averaging the fluxes 
allows for a more accurate representation of the RankineTIugoniot jump conditions in 
one dimension. In this paper we average the fluxes but have not seen much difference 
between the various approaches. 

The dissipation terms are a blending of second and fourth differences. That is, 


AD = (d* + D 2 v - D\ - D*) W, 


(5) 


where 


D \ W = A <] 


D\W = V i 


, A4) 




w 

1>3 ? 


( 6 ) 

(7) 


and A^ , are the standard forward and backward difference operators, respectively, 
associated with the £ direction. The variable scaling factor A is chosen as 


K U+l 


1 

2 




+1J 


( 8 ) 


where A^ is proportional to the spectral radius of the matrix A. The coefficients e< 2) and 
are adapted to the flow and are defined as follows: 


e (2) 

*+W 


= *( 2 ) 


max(i/,. lt j, Vij, I'i+ij, v t + 2j), 


(9) 
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( 10 ) 
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Pi+lJ - ‘^PPj ± Pi-l ,j 
Pi+ 1 ,j T 2p,-,j T Pi—ij 


max 


0, 



( 11 ) 


where p is the pressure, and the quantities and are constants to be specified. 
The operators for the p direction are defined in a similar manner. 

The second-difference dissipation term is nonlinear. Its purpose is to introduce an 
entropy-like condition and to suppress oscillations in the neighborhood of shocks. This 
term is small in the smooth portion of the flow field. The fourth-difference dissipation 
term is basically linear and is included to damp high-frequency modes and allow the 
scheme to approach a steady state. Only this term affects the linear stability of the 
scheme. Near shocks it is reduced to zero. For high-speed flows, the switch (10) is not 
very good and does not allow the multigrid to converge. Instead we considei a T\ D 
variation of the switch [11] given by 


Vi,i = 


|pi+l,j ~ + P«-l,jl 

|pi+i,j - Pi.il + I Pt,j - Pt-I.il + e ’ 


= 1 / 2 . 


( 12 ) 


With this change and the factor 1/2 in front of the second-difference dissipation term, 
the scalar equation becomes first-order upwind near shocks. In the case of the original 
{/, we find that v ~ 0.05 near shock waves in transonic flows. The parameter e must be 
chosen carefully to prevent the switch from being activated by noise. In fact, we found 
it useful to take an average of the two versions for v. Hence, we use 


|P»+i.j 2pi,j 4~ Pt-l,j I 

(1 — U?) ( P TV D ) «, j + uPi 


(13) 


where 

( Ptvd)\,j — |Pi+i,i — Pi.jl + IPi.i — Pi-i, il> 

P\,j = Pi+i,i T 2p,,i T Pi— l ,i •> 

and u is taken to be 1/2. We now no longer have a free parameter for the second- 
difference dissipation. 

Several other changes were made to the scheme in addition to the change to a TVD 
switch. In the original algorithm, the artificial viscosity for the energy equation was 
based on the total enthalpy rather than the total internal energy. For high-speed flows 
we base the artificial viscosity on the total internal energy so that in each equation the 
basic dependent variable is also used in the artificial viscosity. This is more in line with 
upwind schemes. This has previously been used in central-difference schemes [14]. The 
algorithm no longer preserves a constant total enthalpy in the steady state (as the Euler 
equations do), but enthalpy damping is not useful for supersonic flows. In most cases the 
difference between the two approaches is small with each approach having its advantages. 
The original version seems to give slightly sharper shocks, while the other one appears 
to make the scheme more robust. 

The form of the dissipation model of the basic or driving scheme is usually modified 
for the coarse grids in the multigrid process. A constant coefficient, second-difference dis- 
sipation is not only less expensive computationally but also generally provides adequate 
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smoothing properties. For high-speed flows we find it necessary to append a nonlinear 
dissipation to the usual one that depends on the modified switching function of (12). We 
also need to increase the constant coefficient on a coarse mesh from the standard value 
of 1/16 to a value of 1/4. 

In order for the scheme to be stable it is necessary to restrict the time step. Both a 
convection limit as well as a diffusion limit must be taken into account in general. As 
shown in [15] the actual time step (Af act ), based on a sufficient condition for stability, is 
determined as follows: 

Af ac t ^ N + A,, + dX viscous ] , (14) 

where N is taken to be the allowable CFL number, and the constant d is 4. The spectral 
radii A^ and A^ of the Jacobian matrices A and 5, where the tilde denotes that the 
matrix is multiplied by the transformation Jacobian J, are given by 

= l u £r + v£ y \ + Cyjtl + 

X v = | VT] y + UT] X \ + CyJrjl + 7 / 2 . 

For the thin-layer Navier-Stokes equations, 


A viscous 


Prp 


(Til + Vi), 


where 7 is the specific heat ratio, and Pr is the Prandtl number. 

We now consider how the contributions to (14) behave near the body as the inflow 
Mach number is increased. We assume that the variables are normalized so that the 
density and pressure are 1 at inflow. In these nondimensional units the viscous terms 
are multiplied by °° . In the boundary layer the velocity goes to zero and so A^ 
basically behaves as the sound speed c. If p = T ~ p/ p, then 


A viscous 

A rj 


Mooy/p 

3 

piRe^Pr 


We next consider the ratios pt 2 /pi and Pt 2 / Pi for adiabatic flow [16]. As M becomes 
larger p behaves like 6 M ^/5 while p approaches 6 . So as M ^ becomes larger 


A viscous j, *2 

— — ~ lVl 

\ ^ J oo 

Apj 

Thus, the diffusion limit on the time step is quite significant for high Mach numbers. 

For all flow calculations in this paper a five stage Runge-Kutta scheme with a weighted 
evaluation, as detailed in [17], of the dissipation terms on the first, third, and fifth stages 
is used. The time step is reduced near shocks by including a term that depends on v^ y 
The reduction is constructed so that there is a CFL number of 1 when v = 1 . It serves to 
reduce the magnitude of the change in the solution near the shock wave, which exhibits 
strong nonlinear behavior. 
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Implicit Residual Smoothing 

A mathematical step can be included in the basic time-stepping scheme to extend 
the local stability. This procedure is called implicit smoothing (or averaging) of the 
residuals. Previously, variable coefficients for this implicit process have been introduced 
by Martinelli [2] and Swanson and Turkel [18]. In addition, an explicit multidimensional 
derivation of such coefficients is considered in [19]. These variable coefficients have 
proven to be quite reliable in extending the stability limit of the Runge-Kutta scheme by 
a factor of two. However, their development is based on hyperbolic considerations only, 
and they do not eliminate the need for a diffusion factor in the time step. As indicated 
in the preceding section and we emphasize here, this factor can even dominate the time 
step determination in the case of hypersonic flow. In this section we will briefly discuss 
the smoothing coefficients of [18] and present their straightforward extension to three 
dimensions. Then we will develop smoothing coefficients that allow the removal of the 
diffusion limit, resulting in a At which depends only on the spectral radii of the flux 
Jacobian matrices. 

For multidimensional problems, the residual smoothing can be applied in the form 


ii(/-av,ai) 


K. 


(m) 


( to ) 


TV 

,v *.j 


/ = £,*?, c 


(15) 


where the residual 7 is defined by 


- a m ^[fc^r 11 + C D W^ - AD^], 


(16) 


and computed in the Runge-Kutta stage m. The quantity AD^ m ^ is the total artificial 
dissipation at stage m, and is the final residual at stage m after the sequence of 

smoothings involving each coordinate direction. The difference operators Cc and Cd 
are associated with the convection and physical diffusion terms. To derive the maximum 
stability extension for the hyperbolic problem, the implicit procedure is applied after 
each stage of the Runge-Kutta scheme. The coefficients are functions of the spectral 
radii A^, A n , and A^. In two dimensions they are written as follows: 


ft 


A? 


max < — 
4 


max < — 

I 4 


N 


1 


k N* 1 + i>r r)( . i 


N 


1 


, N * 1 + 



(17) 

(18) 


where the ratio r v { = A^/A^, and the quantity N/N* is the ratio of the CFL number 
of the smoothed scheme to that of the basic explicit scheme (usually having a value of 
2). In hypersonic flow applications we found it necessary for N* to be 3.25, rather than 
the value of 3.75 used for transonic computations. From a linear stability analysis, the 
scheme with these coefficients is stable for all mesh cell aspect ratios when the parameter 
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ip ~ 0.125 and N/N* is sufficiently large. The practical limitation on the Courant number 
is due to the requirement for effective high frequency damping. For large N/N*, the 
high frequency damping of the scheme vanishes. The variable coefficients are functions 
of the local mesh cell aspect ratio, and thus the smoothing process is not activated 
in a coordinate direction where it is not needed. This is important for best possible 
convergence. 

The formulas of (17) and (18) can be easily extended to three dimensions. Moreover, 
we define 


where 


Pi = max l - 


N ' 2 

— 

N* 


r- 


.0 




1 + V>3-r>(»V + r cc )’ 


1 + ipz -D(r n l + r^)’ 


1 + V*3 -o( r c/ + r (t] ) 


(19) 


Once again r represents a ratio of characteristic speeds. A typical value for is 

0.0625. 

We now show how one can utilize these coefficients as the basis for new coefficients 
that will remove the diffusion limit on the time step. Consider the scalar diffusion 
equation 

dw 8 2 w 

lit = 

If we approximate the spatial derivative of (20) with a central difference and take a 
Fourier transform, we obtain 

At^-=Zw n , (21) 

where the caret indicates a transformed quantity, and the Fourier symbol 



Z — — 2Ad(l — cosd ) 


( 22 ) 


with A d = Atfj,/ Ay 2 . If implicit residual smoothing is applied, then 


Z 


— 2Aj(l — cos0) 

1 +2/? d (l -co50) 

A sufficient condition for stability is given by 

|Z| < N* d for all 0, 


(23) 


(24) 


where N^ is the diffusion number of the unsmoothed scheme. It then follows that the 
smoothing coefficient (3 d is given by 


Pd>\ 

A 


L^f . 


A u 

At* 


- 1 


(25) 
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where Nd — 4A This implicit smoothing coefficient has a linear dependence on the dif- 
fusion number ratio. In the case of the pure hyperbolic problem, one obtains a quadratic 
dependence with the CFL number ratio. Thus we must determine a way to combine 
these results to obtain a coefficient that is valid for an equation with both convection 
and diffusion effects. 

Suppose we consider the thin-layer form of the two-dimensional Navier-Stokes equa- 
tions, then one could simply use the smoothing coefficient of (17) in the streamwise-like 
(£) direction. One would anticipate that a possible formulation for the normal r/ direc- 
tion would depend on a diffusion type (3 near the surface and a convection type (3 when 
the viscous effects are no longer important. Notice the dependency in (25) on the ratio 
of the actual At to the At of the basic explicit scheme. So we define 


w. 


At 


act 


[(A^)r? 


- 1 


(26) 


Since we want to remove the diffusion limit, the actual At must depend only on A^ and 
A ?r Thus we rewrite (26) as 




1 

4 


Cx 


A^ + A, 


-1 , 


which we can approximate simply as 




__ m (Mt; 

4 1 + A/ 


(27) 


(28) 


It does not appear to matter whether one uses the form of (27) or that of (28), provided 
the constant is defined properly. Here we have elected to use the simpler form of (28), 
which has also been considered by Radespiel and Kroll [20]. From numerical experiments 
we found that a satisfactory value for C\ is 10. 

The variable coefficient of (28) cannot generally be used by itself. For example, in 
an airfoil flow, (f3d) v goes to zero too fast at the leading edge, resulting in a zero value 
in the inviscid region. We have overcome this difficulty by calculating (3 V as 


/?„ = max ((/?<*)„ , (0 C ) u) , (29) 

where (0 c )j? is defined by (18). In the Results section of this paper we consider Mach 10 
and Mach 20 turbulent flow over an airfoil. The removal of the diffusion limit allowed 
the residual in each case to be reduced at least an additional order of magnitude in 200 
multigrid cycles. 


Multigrid Method 

As indicated earlier the salient features of the multigrid method considered here are 
fairly standard. We apply the Full Approximation Storage (FAS) scheme of Brandt [21] 
to define the equivalent fine grid problem on a coarse grid. The grid transfer operators 
for the solution, residual, and coarse grid corrections are those introduced by Jameson 
[1], In order to execute the multigrid strategy we employ a fixed W-type cycle with two 
sweeps on each coarse grid. To provide a well conditioned starting solution for the fine 
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mesh a Full Multigrid (FMG) method is used. The FMG is analogous to grid sequencing, 
except that multigrid cycles are performed on each coarse grid. 

Some of the additional elements of the multigrid method are not necessarily standard. 
A smoothing of the coarse grid corrections being transferred to the finest grid was found 
to be beneficial in transonic computations [17]. The smoothing was accomplished with 
the implicit residual smoothing mentioned previously and a constant coefficient /? ~ 0.1. 
This smoothing of the residuals on the way to finer meshes is crucial for the convergence 
of the multigrid for hypersonic flows. Such a process acts to reduce high frequency oscil- 
lations caused by the interpolation. Hence, it becomes especially important near strong 
shocks, where nonphysical upstream influence can occur. It should be emphasized that 
choosing the smoothing parameter too large can slow convergence. Another important 
element for high Mach number (M >10) flows is the coarse grid correction scheme. That 
is, the physical viscous terms should also be computed on the coarse meshes. Difficulties 
with applying a turbulence model on a coarse grid can be avoided by interpolating the 
turbulent viscosity from the values on the finest grid. 

Boundary Conditions and Initialization 

At a solid surface (wall) boundary the no-slip condition is enforced. The wall pres- 
sure is set to the value at the first interior solution point, and thus, a reduced normal 
momentum equation is satisfied. An adiabatic wall is assumed. In a finite-volume for- 
mulation, this amounts to treating the Cartesian velocity components as antisymmetric 
functions and the temperature as a symmetric function with respect to the wall. For each 
of the physical problems considered the Mach number at the inflow boundary exceeds 
1.0. Consequently, the dependent variables are specified at this boundary according to 
the flow conditions. At any outflow boundary, we apply simple extrapolation of the 
components of the solution vector. In general, for hypersonic flows, numerical difficulties 
are experienced at the start of a calculation if the discrete flow field is initialized with 
free-stream conditions. To avoid these difficulties we apply the following procedure. The 
Mach number of the flow is set to a lower value than the required one. Then the Mach 
number is gradually increased over a few hundred time steps until the desired flow con- 
ditions are obtained. This Mach number ramping is only done on the coarsest mesh in 
the FMG sequence. 

Results 

We consider the following two test cases: (1) two-dimensional turbulent flow about 
the NACA 0012 airfoil, and (2) three-dimensional turbulent flow about a blunt biconic. 
For these cases we found a significant sensitivity to the grid. With a grid suitable for 
transonic flows (i.e., a standard C-type mesh), difficulties in convergence were expe- 
rienced at high Mach numbers. Furthermore, the usual mesh density for airfoil flows 
yielded a poor resolution of the flow field. The grids which were generated for these 
examples all follow in a general sense the bow shock. The bow shock is not aligned with 
the grid, but as a minimum, the angle between the bow shock and the grid is not too 
large. We also found it important to control the normal stretching rate of the grid in 
the inviscid region of the flow field. 

For the first case we consider turbulent flow (a Reynolds number Re c — 10 7 ) about 
the NACA 0012 airfoil at two different Mach numbers. In all the calculations, transition 


9 


was assumed to occur at the 20 percent chord location. This ensured that the algebraic 
turbulence model of Baldwin and Lomax [22] did not determine an unrealistic length 
scale in the vicinity of the leading edge due to the presence of the shock. The grids for 
the computations were generated with the method of [23] . The fine grid, consisting of 
320 x 64 cells, is displayed in figure 1. In the streamwise direction, the grid is clustered at 
the leading and trailing edges of the airfoil, with a tangential spacing of approximately 
1 X 10 -3 chords. There is clustering in the normal direction adjacent to the surface, so as 
to provide adequate resolution of the turbulent boundary layer. The minimum normal 
spacing occurs at the airfoil leading edge and is 1 x 10~ 5 chords. To obtain reasonable 
shock definition everywhere, there is weak normal clustering in the inviscid region. 

Figures 2 and 3 show the Mach number and pressure contours, respectively, when 
Moo = 10. There is fairly good resolution of the bow shock, and the contours are 
smooth. The pressure contours delineate the expected fish tail shock in the wake region. 
Surface pressure and skin-friction distributions for this case are presented in figure 4. 
In addition, velocity profiles at the midchord and 95 percent chord locations are shown 
in figure 5. The same results computed on a 160 x 32 mesh, which is a proper subset 
of the 320 x 64 mesh, are also included in these figures. One can clearly see that there 
is fairly good agreement between the predictions on the two meshes. Thus, additional 
mesh refinement will serve primarily to improve the downstream resolution of the bow 
shock. In figure 6 the variations of the Mach number and pressure along the stagnation 
streamline are shown. Comparing the pressure jump and the stagnation pressure with 
one-dimensional theory, we find the differences are less than 2 percent. With the cell- 
centered, finite- volume scheme, one should keep in mind that the computational points 
do not lie precisely along the stagnation streamline. 

Convergence histories for the Mach 10 case are given in figure 7. On the finest 
mesh the logarithm of the residual of the continuity equation is reduced 7.5 orders of 
magnitude in 200 multigrid cycles. This corresponds to an average rate of reduction of 
0.917. The calculation required less than 6 minutes on a Cray II computer. It should 
be emphasized that for engineering accuracy (i.e., residual reduced by 3 orders) the 
finest mesh calculation required less than 3 minutes of CPU time. Note that when 
engineering accuracy is achieved, there is no appreciable improvement in the viscous 
solution accuracy by further residual reduction. 

In figures 8-11 results for = 20 and Re c = 10 7 are presented. Again there is fairly 
good representation of the bow shock, and the Mach and pressure contours are smooth. 
Here also there is less than a 2 percent difference between the computed pressure jump 
and the stagnation pressure determined along the stagnation streamline and the values 
from one-dimensional theory. As indicated in figure 12, the logarithm of the residual 
is reduced nearly 5 orders in 200 cycles on the finest mesh. There is a slowdown in 
the asymptotic convergence rate for this case. At this point the precise reason for such 
behavior is not clear. 

We now consider the three-dimensional case. A mesh consisting of 128 x 96 x 24 
cells was used for the blunt biconic. In figure 13, a two-dimensional slice of the mesh is 
depicted. This was constructed using an algebraic mesh generator. A tangent hyperbolic 
distribution function was used in the normal direction to obtain an acceptable mesh in 
both the viscous and inviscid regions. For this turbulent flow case, M^ = 6, a — 5 
degrees, and the Reynolds number based on nose diameter is 2.89 x 10 5 . Figures 14 and 
15(a) - 15(b) show Mach number and pressure contours, respectively, in the symmetry 
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plane. One sees the small standoff distance of the shock at the nose and how the shock 
remains relatively close to the geometry downstream. The shock surface is delineated 
in the three-dimensional view of pressure contours given in figure 15(c). A comparison 
is made in figure 16 between computed pressure distributions, for the windward and 
leeward planes, and the experimental data of [24]. There is generally good agreement 
with the data. In figure 17 the convergence history is presented for this computation. 
The logarithm of the residual has been reduced about 5.5 orders of magnitude in about 
565 work units, which roughly corresponds to 350 multigrid cycles on the finest mesh. 
For engineering accuracy, 85 cycles are required on the finest mesh, and the CPU time 
is about 30 minutes. 

Concluding Remarks 

A multigrid method with central differencing has been successfully applied to the 
solution of hypersonic viscous flows. An explicit five stage Runge-Kutta scheme has been 
used as a smoother in solving the time-dependent, thin-layer Navier-Stokes equations. 
In this paper, considerable emphasis has been focussed on the dissipative characteristics 
of the driving scheme for the multigrid process. The presence of strong shocks has 
required the introduction of a switching function for the numerical dissipation based 
on TVD principles. This nonlinear switching function is also required for the coarse 
grid dissipation model. The importance of the physical diffusion limit on the time step 
for computing hypersonic flows has been discussed. By introducing appropriate variable 
coefficients for implicit residual smoothing, we have removed the diffusion limit and have 
improved the convergence significantly. 

Numerical solutions have been obtained for two- and three-dimensional hypersonic 
turbulent flows. The agreement between predictions for flow over a blunt biconic compare 
well with experimental data. Engineering accuracy has been obtained rapidly in all 
computations, requiring less than 3 CPU minutes on a Cray II for two-dimensional cases 
and about 30 CPU minutes for a three-dimensional case. 
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(a) 320 x 64 grid 


(b) Blowup of leading edge region 


Figure 1 A grid for turbulent flow over NACA 0012 airfoil 



Figure 2 Mach number contours for Mach 10 turbulent 
flow over NACA 0012 airfoil (a = Odeg, Re c = 10 7 ) 
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Figure 5 Velocity profiles for Mach 10 turbulent 
flow over NACA 0012 airfoil (a = Odeg, Re c = 10 7 ) 
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Figure 6 Pressure and Mach number variations 
along stagnation streamline for Mach 10 turbulent 
flow (NACA 0012 airfoil, a = 0 deg, Rt c — 10 7 ) 


Figure 7 Convergence histories for Mach 10 turbulent 
flow (NACA 0012 airfoil, a = 0 deg, Re c = 10 7 ) 
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(a) Surface pressure coefficients 


(b) Surface skin-friction coefficients 


Figure 10 Surface pressure and skin-friction coefficients for Mach 
20 turbulent flow over NACA 0012 airfoil (a = 0 deg, Re c = 10 7 ) 
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Figure 11 Pressure and Mach number variations 
along stagnation streamline for Mach 20 turbulent 
flow (NACA 0012 airfoil, a = 0 deg, Re c = 10 7 ) 


Figure 12 Convergence histories for Mach 20 turbulent 
flow (NACA 0012 airfoil, a = 0 deg, Re c — 10 7 ) 
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(a) Windward plane 


(b) Leeward plane 


Figure 16 Pressure coefficient distributions in axial direction for Mach 
6 turbulent flow over blunt biconic (o = 5 deg, Rep - 2.89 x 10 5 ) 



Figure 17 Convergence history (blunt biconic, a = 5 deg, Re p - 2.89 x 10 5 ) 
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